The pressure distribution in thermally bistable turbulent flows 
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ABSTRACT 

We present a systematic numerical study of the effect of turbulent velocity 
fluctuations on the thermal pressure distribution in thermally bistable flows. The 
turbulent fluctuations are characterized by their rms Mach number M (with re- 
spect to the warm medium) and the energy injection wavenumber, kf or = l/£, 
where I is the injection size scale in units of the box size L =100 pc. The nu- 
merical simulations employ random turbulent driving generated in Fourier space 
rather than star-like heating, in order to allow for precise control of the parame- 
ters. Our range of parameters is 0.5 < M < 1.25 and 2 < kf or < 16. Our results 
are consistent with the picture that as either of these parameters is increased, the 
local ratio of turbulent crossing time to cooling time decreases, causing transient 
structures in which the effective behavior is intermediate between the thermal- 
equilibrium and adiabatic regimes. As a result, the effective polytropic exponent 
7 e of the simulations ranges between ~ 0.2 to ~ 1.1, and the mean pressure of 
the diffuse gas is generally reduced below the thermal equilibrium pressure P eq , 
while that of the dense gas is increased with respect to P oq . The fraction of 
high-density zones (n > 7.1 cm -3 ) with P > 10 4 K cm -3 increases from roughly 
0.1% at k {or = 2 and M = 0.5 to roughly 70% for k {or = 16 and M = 1.25. A pre- 
liminary comparison with the recent pressure measurements of Jenkins (2004) 
in CI favors our case with M = 0.5 and /cf or = 2. In all cases, the dynamic 
range of the pressure in any given density interval is larger than one order of 
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magnitude, and the total dynamic range, summed over the entire density range, 
typically spans 3-4 orders of magnitude. The total pressure histogram widens 
as the Mach number is increased, and moreover develops near-power-law tails at 
high (resp. low) pressures when 7 e < 0.5 (resp. 7 C > 1), which occurs at kf or = 2 
(resp. fcfor = 16) in our simulations. The opposite side of the pressure histogram 
decays rapidly, in an approximately lognormal form. This behavior resembles 
that of the corresponding density histograms, in spite of the large scatter of the 
pressure in any given density interval. Our results show that turbulent advec- 
tion alone can generate large pressure scatters, with power-law high-P tails for 
large-scale driving, and provide validation for approaches attempting to derive 
the shape of the pressure histogram through a change of variable from the known 
form of the density histogram, such as that performed by Mac Low et al. (2004). 

Subject headings: ISM: structure — instabilities — turbulence — hydrodynamics 
— ISM: kinematics and dynamics 

1. Introduction 

The atomic interstellar medium (ISM) is generally believed to be thermally bistable. 
This property arises because the neutral gas is thermally unstable for 300 K < T < 5000 K 
under the isobaric mode of thermal instability (TI; Field 1965; see also the review by Meerson 
1996), allowing for a configuration in which gas with temperatures above and below this range 
can coexist in thermal pressure equilibrium (Field, Goldsmith & Habing 1969; Wolfire et al. 
1995, 2003), being mediated by a thin interface of thickness comparable to the conductive 
length (Begelman & McKee 1990). This tendency of the HI gas to naturally segregate in two 
phases has long been thought to be the dominant mechanism in forming and maintaining 
cold cloudlets of sizes ~ 0.1 pc confined by the thermal pressure of a warm, dilute substrate. 

On the other hand, the ISM is known to be globally turbulent, with numerous kinds of 
energy sources acting on a wide range of scales (e.g., Scalo 1987; Norman & Ferrara 1996), 
including spiral arm shocks, large-scale combined instabilities, supernova and nil-region 
energy input, etc. (see, e.g., Mac Low & Klessen 2004 for a review). The net effect is to 
produce and maintain a turbulent velocity dispersion that is transonic with respect to the 
warm gas, and supersonic with respect to the cold medium (e.g., Heiles & Troland 2003). 

In recent years the interplay between TI and turbulence has been studied by several 
groups. Some of the main issues are: The development and signatures of TI in a turbulent 
ISM (e.g. Vazquez-Semadeni, Gazol & Scalo 2000; Sanchez-Salcedo, Vazquez-Semadeni & 
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Gazol 2002; Vazquez- Semadeni et al. 2003), the triggering of TI by strong compressions 
in the warm medium (Hennebelle & Perault 1999, 2000), the generation of turbulence at 
the nonlinear stages of development of TI (Koyama & Inutsuka 2002; Kritsuk & Norman 
2002), the production and maintenance of small-scale cold structures (Burkert & Lin 2000; 
Koyama & Inutsuka 2000, 2002; Audit & Hennebelle 2005), the interaction between TI and 
magneto-rotational instability (Piontek & Ostriker 2004) the production and maintenance 
of significant amouts of gas at thermally unstable temperatures (Gazol et al. 2001; Audit 
& Hennebelle 2004), as motivated by the reports of several observational studies (Dickey, 
Salpeter & Terzian 1977; Kalberla, Schwarz & Goss 1985; Spitzer & Fitzpatrick 1995; Fitz- 
patrick & Spitzer 1997; Heiles 2001; Heiles & Troland 2003; Kanekar et al. 2003), and the 
thermal pressure distribution in the turbulent ISM (Audit & Hennebelle 2004; de Avillez & 
Breitschwerdt 2004; Mac Low et al. 2004; see also Vazquez-Semadeni et al. 1995; Korpi et 
al. 1999 for discussions of cases without thermal bistability.) 

Physical discussions of the effect of turbulent velocity fluctuations in a thermally bistable 
medium have been given by Sanchez- Salcedo & al. (2002), Vazquez-Semadeni et al. (2003), 
Wolfire et al. (2003), and Audit & Hennebelle (2004). The first two of these works noted 
that velocity fluctuations induce perturbations on the gas that can range from behaving 
adiabatically, when the turbulent crossing time r t across the perturbation size scale is much 
shorter than the cooling time r c , to behaving according to the thermal equilibrium condition 
between cooling and heating, in the opposite limit. The turbulent crossing time, in turn, 
depends on the scale and amplitude of the velocity fluctuations, and therefore the higher the 
Mach number, or the smaller the typical scale of the turbulence, the higher the fraction of 
fluid parcels that are expected to transiently behave closer to an adiabatic regime, and farther 
away from thermal equilibrium. Sanchez-Salcedo & al. (2002) and Vazquez-Semadeni et al. 
(2003) used this to explain the presence of significant amounts of gas with temperatures 
corresponding to the unstable range. Audit & Hennebelle (2004) have further quantified the 
problem by producing a semi-analytical model that follows the stretching of a fluid parcel 
due to the shearing components of the turbulence while it seeks thermal equilibrium, to give 
a relation between the amplitude of this component and the amount of thermally unstable 
gas present in the flow. Finally, Wolfire et al. (2003) have given an estimate of the ratio rj of 
the turbulent crossing time to the cooling time in the warm neutral medium, finding values 
0.3-0.9, which led them to suggest that this medium should often exhibit non-equilibrium 
temperatures. 

These results have implications for the thermal pressure PDF in the flow. As a fluid 
parcel departs from thermal equilibrium, its pressure also departs from the equilibrium value, 
and we expect a distribution of the thermal pressure around its thermal equilibrium value at 
a given density, determined by the distribution of Mach numbers of the velocity fluctuations. 
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This holds even in the absence of direct local heating. 

The thermal pressure probability density function (PDF) varies significantly among 
different models of the ISM. In the simplest equilibrium multiphase model (Field et al. 1969), 
the thermal pressure PDF is simply a delta function at the mean value of the pressure in the 
midplane. If this value encompasses the unstable range, then the medium is expected to be 
segregated into the two phases, both at the mean, equilibrium pressure, but with densities 
and temperatures bracketing the unstable range. The next level of complexity was added by 
including supernova heating (Cox & Smith 1974; McKee & Ostriker 1977). In particular, the 
McKee & Ostriker (1977) model implied a piecewise power-law pressure PDF (Jenkins, Jura 
& Lowenstein 1983; see also Mac Low et al. 2004) with slope -19/9 for P < P c and -23/9 
for P > P c , where P c /k = 10 3 67 cm~ 3 i^. Moreover, the model predicted no pressures below 
the equilibrium pressure of the warm and cold phases. This distribution, however, follows 
from a consideration of the probability of a given point in space belonging to a hot supernova 
remnant, and the theoretical evolution of these remnants, not from a consideration of the 
local thermodynamic changes in the gas due to the turbulent compressions and rarefactions 
(advection) that are induced by the supernova energy injection. 

Inclusion of advection is naturally accomplished in numerical models of the star-driven 
ISM (see Vazquez-Semadeni 2002 for a review). In particular, the recent papers by Mac 
Low et al. (2004) and de Avillez & Breitschwerdt (2004) have discussed the pressure PDF 
resulting in their simulations, albeit they appear to obtain different functional forms for it: 
Mac Low et al. find a lognormal PDF, while de Avillez & Breitschwerdt find PDFs that 
appear closer to a power law, with a slope in fact not too different from that predicted 
by McKee & Ostriker (1977). From the lognormal shape of their PDF, Mac Low et al. 
conclude that it originates from the density PDF for an isothermal turbulent flow. However, 
in those simulations it is not possible to disentangle the pressure fluctuations induced purely 
by turbulent motions, and those due to direct heating from nearby stellar sources. 

Observationally, significant pressure fluctuations have been reported in the cold medium. 
Jenkins et al. (1983) found, using Copernicus observations of CI, variations greater than an 
order of magnitude in the cold gas pressure with small amounts of gas at up to P/k = 
10 5 K cm -3 . More recently, Jenkins & Tripp (2001) used the Space Telescope Imaging 
Spectrograph (STIS) to confirm this result with better resolved data. Additionally they 
found that their results implied an effective polytropic index for the cold gas 7 > 0.9, which 
is larger than the 7 = 0.72 derived by Wolfire et al. (1995) for cold gas at thermal equilibrium, 
mentioning that this could be due to the fact that compressed regions may have a cooling 
time larger than the dynamical time, so it may behave closer to adiabatically. Finally, it is 
well known that in the local ISM there is an apparent pressure imbalance between the hot 
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(T ~ 10 6 K, P/k b ~ 11, OOOK cm" 3 ) and the warm (T ~ 6, 700K, P/k b ~ 2, 300K cm- 3 ) gas 
(see e.g. Jenkins 2002; Redfield & Linsky 2004). 

In this paper, we present a systematic study aimed at investigating in detail the ef- 
fects of turbulent velocity fluctuations on a thermally bistable flow, in particular on the 
transition from nearly-adiabatic to near-thermal equilibrium behavior, and on the PDF of 
thermal pressure. The simulations we present here are not intended as accurate models of 
the ISM, but instead as numerical experiments allowing us to clearly identify the effects of 
two fundamental turbulent parameters, the rms Mach number M and the energy injection 
scale (characterized by its wavenumber fcf or ), on the thermodynamic response of a thermally 
bistable flow. For this reason, we have opted for using Fourier random driving, trading 
up realism for accurate control of these parameters, and consider non-magnetic, non-self- 
gravitating flows. Also, in order to allow for the large number of simulations needed to cover 
a significant range in parameter space, we have restricted most of the simulations to two 
dimensions, although we present a few selected cases in three dimensions (3D) in order to 
check that the main trends observed the two-dimensional runs are preserved in 3D. 

The outline of the paper is as follows. In §2 we describe the model used for the simu- 
lations and its limitations, and in §3, we present convergence tests for this model. In §4 we 
then present the main results, concerning the effective thermodynamic behavior of the flow 
and the pressure PDF, both in global form and in specific density intervals, as the parameters 
M (§4.1) and fcf or (§4.2) are varied. Then, in §5 we discuss these results in the context of the 
simple physical scenario of a transition from near-thermal-equilibrium to near-adiabaticity 
(§5.1), we show persistence of the main trends in 3D (§5.2), and we discuss some implications 
of our results for previous models and simulations (§5.3). Finally, in §6 we give a summary 
and some conclusions. 



2. The Model 
2.1. Prescription 

We solve the hydrodynamic equations, including the energy equation, to simulate a 
region of 100 pc on a side, with periodic boundary conditions. The simulations discussed in 
this paper are in two dimensions (2D) except for those described in §5.2. 

The equations are solved by means of a MUSCL-type scheme (Monotone Upstream- 
centered Scheme for Conservation Laws) with HLL Riemann solvers (Harten, Lax, & van 
Leer 1983; Toro 1999), augmented with model terms for radiative cooling and background 
heating, and a prescription for random turbulent forcing. The background heating is taken 
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as a constant T = 2.51 x lCT 26 erg s _1 H _1 , where "H -1 " means "per Hydrogen atom". This 
is the value of the photo-electric heating rate at density n = 1 cm~ 3 reported by Wolfire et 
al. (1995), and is roughly within half an order of magnitude of its value throughout the range 
10~ 2 cm~ 3 < n < 10 3 cm~ 3 , as reported by those authors. We then use this value to fit the 
"standard" equilibrium P vs. p curve of Wolfire et al. (1995) assuming that the background 
heating is in equilibrium with a cooling function that has a piece-wise power-law dependence 
on the temperature. We find (Sanchez-Salcedo et al. 2002) 



A 



T < 15 K 

3.42 x io 16 T 213 15 K < T < 141 K 

9.10xl0 18 T 141 K < T < 313 K 

1.11 x 10 20 T - 565 313 K < T < 6101 K 

2.00 x 10 8 T 3 67 6101 K < T 



where the coefficients have units of erg s~ 1 g~ 2 cm 3 K~ /3 , with (5 being the temperature ex- 
ponent in the corresponding interval. Under thermal equilibrium conditions, the gas is 
thermally unstable under the isobaric mode for 313 K < T < 6102 K, and marginally stable 
for 141 K < T < 313 K. In thermal equilibrium, the transition temperatures T = 6102, 313 
and 141 K correspond to densities n = 0.60, 3.2 and 7.1 cm -3 , respectively. 

The turbulent driving is 100% solenoidal, and is done in Fourier space at a specified 
narrow two- or three-dimensional wavenumber band, kf or — 1 < k < k{ or , where k = ^k"l + ky 
in 2D and k = ^k\ + A; 2 + k 2 z in 3D, and with the Gaussian deviates having zero mean and 
unitary standard deviation. The amplitude of velocity perturbations is fixed by a constant 
injection rate of kinetic energy as in the prescription of Mac Low (1999), although with the 
difference that we use a different random seed at each driving time. The kinetic energy input 
rate is chosen as to approximately maintain a desired sonic Mach number. 

In the set of simulations presented in §4, the fluid is initially at rest and has a uniform 
density (no = 1 cm -3 ) and temperature (T = 2399 K), so that, in the absence of turbulence, 
the medium would spontaneously segregate into warm-diffuse (n = 0.34 cm -3 , T = 7104 K) 
and cold-dense (n = 37.2 cm -3 , T = 64.5 K) phases. The time unit to is chosen to be the 
sound crossing time across the numerical box at a speed of 9.1 km s _1 , corresponding to the 
isothermal sound speed at 10 4 K. Thus, t = 10.8 Myr. Mach numbers are expressed with 
respect to this sound speed. 

For the parameters we use as initial conditions, and in the presence of realistic thermal 
conductivity (see §3), the maximum linear growth rate occurs at scales ~ 8.3 pc, while 
the so called Field length (the minimum unstable scale) is 0.7 pc. In our simulations, 
thermal conductivity is not included, other than the numerical diffusion caused by the finite- 
differencing. This causes a "numerical Field length" ~ 3 pc in simulations at a resolution 
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of 512 2 , with the maximum growth rate occurring at scales between 12.5 and 25 pc, as 
determined through the tests described in §3. Thus, the unstable wavelength range in our 
simulations is somewhat more compressed than the real unstable range in the linear regime. 
We discuss the neglect of thermal conductivity further in §2.2.2, and resolution issues in §3. 

2.2. Features and limitations of the model 

2.2.1. Random Fourier driving 

In the simulations presented here, the turbulence is driven using a random scheme 
executed in Fourier space. This means that every point in physical space is subject to a 
force at any given time. This is not a very realistic way of driving the turbulence, as in 
the real ISM the driving sources, such as supernova explosions or spiral-arm shock waves 
are localized in space. Nevertheless, turbulence is expected to "propagate away" from the 
localized sources (Avila-Reese & Vazquez- Semadeni 2001), generating a general turbulent 
flow. More importantly, we choose this form of driving because of two reasons. First, it 
allows us to precisely control the scale of energy injection as well as the rms Mach number 
in the flow. This is very important, because one of our motivations is to test the scenario 
described in Sanchez-Salcedo et al. (2002) and Vazquez-Semadeni et al. (2003) that velocity 
fluctuations behave closer to an adiabatic regime as the turbulent crossing time becomes 
shorter than the cooling time, and the turbulent crossing time is directly a function of the 
scale size and amplitude of the velocity perturbations. 

The second reason is that random Fourier driving guarantees that all pressure fluc- 
tuations that develop are caused purely by advection (fluid transport) and not by direct 
injection of heat by stellar sources, allowing us to isolate the effects of velocity fluctuations 
on the pressure distribution. This implies that the widths of the pressure distributions in our 
simulations should constitute a lower limit to the actual widths expected in the actual ISM, 
in which the energy is injected as heat by stellar sources, directly raising the local pressure, 
in addition to any effects of the turbulent advection. 

2.2.2. Neglect of thermal conductivity 

The equations solved in our simulations do not include a model term for the thermal 
conductivity. Koyama & Inutsuka (2004) have suggested that thermal conductivity should 
always be included in numerical simulations of thermally unstable flows, and that enough 
resolution to resolve the so called "Field length" (Field 1965; Begelman & McKee 1990) 
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should be always be used, because of mainly two reasons. First, if no thermal conduction at 
all is included, then the TI growth rate asymptotically approaches its maximum value as the 
perturbation length scale approaches zero. This means that, in a finite-resolution numerical 
grid, the smallest resolved scale is maximally unstable, regardless of the resolution used, 
creating numerical problems. However, in practice this does not occur, as the numerical 
diffusion of the code creates a "numerical Field length" , that is, a minimum unstable scale 
larger than the grid cell size, even if it does not have the same functional temperature 
dependence as the actual Field length. Indeed, we have verified that density perturbations 
of amplitude 2.5% and wavelength A = 16 pixels (3.1 pc) in a simulation of resolution 512 2 
remain stationary (thus being the "numerical Field length"), while pertubations with A = 4 
pixels (0.8 pc) are completely damped in times ~ 3 Myr. This timescale is comparable to the 
e- folding time of the growing modes described in §3. Thus, numerical diffusion adequately 
prevents instability of the smallest resolved scales. Because of the very small wavelength, 
numerical dissipation is more effective than TI, which results in damping rather than growing 
of the perturbation. 

Note that Koyama & Inutsuka (2004) also warned that if no realistic thermal conductiv- 
ity is employed, then the results are resolution-dependent, because the characteristic scale of 
numerical diffusion depends on the resolution. Thus, they concluded that, in order to have 
numerically converged results in the presence of TI, one should i) include an explicit conduc- 
tion term that is larger than the conduction due to numerical diffusion, ii) use a cell size that 
is smaller than one third of the conductive Field length. However, for the global properties 
that interest us here, such as PDFs of thermal pressure and density, the difference between 
the effects of numerical diffusion and of an explicit conduction term, whose associated Field 
length is comparable or smaller than the cell size, is probably not significant. Indeed, in §3 
we show that our simulations are perfectly converged at high densities and pressures (which 
is normally the regime of most concern; e.g., Audit & Hennebelle 2004) at the resolution we 
use. 

The second warning of Koyama & Inutsuka (2004) refers to the possibility of missing 
dynamical effects originating from the effects of the thermal conduction. Specifically, they 
describe the generation of motions with Mach numbers up to 0.13 due to the pressure gra- 
dients generated by the conductivity. Moreover, in their simulations of the development of 
TI alone, these motions cause the resulting condensations to coalesce, so that at the end of 
the simulation the number of condensations has decreased by almost a factor of 2. However, 
Koyama & Inutsuka (2004) noted that the initial number of condensations formed in their 
simulations was determined by the initial fluctuations, not by the inclusion or omission of 
thermal conductivity. In addition, the Mach numbers generated in their simulations are at 
least 5 times smaller than those of the smallest turbulent motions we impose on the flow, and 
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therefore, in our simulations, density fluctuation production and coalescence is dominated 
by the turbulent velocity fluctuations, not by TI in the presence of thermal conductivity. 
Thus, for our purposes, the motions produced by thermal conductivity can be neglected. 

3. Convergence test 

Here we present results from 20 simulations tailored to investigate the convergence of 
the model described in §2. 

We first consider the linear regime. In figure 1 we show the temporal growth of sinu- 
soidal density perturbations with an initial amplitude of 2.5%. The spatial period of these 
perturbations is 50 pc (upper left panel), 25 pc (upper right panel), 12.5 pc (lower left panel) 
and 6.25 pc (lower right panel), respectively corresponding to wave-numbers (k p ) of 2, 4, 
8 and 16 respectively. The solid, dotted and dashed lines are for resolutions of N = 256 2 , 
512 2 and 1024 2 , respectively. The thick line in each frame indicates the theoretical growth 
rate at the corresponding scale. To obtain this slope, we have solved the dispersion relation 
(Field 1965) using the cooling function described above (eq. (1)) and a realistic, although 
temperature-independent, value of the conductivity K = K = 5/3(kbT l/v rms )no(3kb/2m) 
(Lang 1999), where k b is Boltzmann's constant, and we have taken T = 2400K, n — 1 
cm -3 , m = m H , I = 3.2 x 10~ 3 pc and v rms = 5.7 km s" 1 , the adiabatic sound speed at T Q . 
These simulations are initially at rest and have a uniform temperature of T = 2400K. It 
can be seen that for k p = 2, 4, and 8 the linear growth rate of the perturbations in the 
simulations at all three resolutions is in good agreement with the theoretical growth rate. 
For k p = 16, we see that the run with N = 512 has nearly converged to the correct growth 
rate, while the run with N = 256 severely damps the growth of this mode. We conclude 
that a resolution iV = 512 is an acceptable resolution for capturing the linear growth, in the 
presence of realistic conductivity, of modes with sizes down to 1/16 the box size. 

We now turn to the nonlinear case, which is the most relevant one for the driven- 
turbulence simulations presented in this paper. In figure 2 we display the temporal growth 
of sinusoidal density perturbations with an initial amplitude of 2.5% for simulations with 
sinusoidal large-amplitude (Mach number 1.0 with respect to the unperturbed medium) ve- 
locity perturbations. The velocity and density perturbations are in phase, with wavelengths 
6.25 pc (left) and 12.5 pc (right). In this case the comparison with the theoretical growth 
rate is not meaningful. However, it can be seen that the difference between the growth rates 
for N = 512 and N = 1024 at k p = 16 is smaller than for the pure density perturbation case. 

As a final test, we compare the total time-averaged (from 1.1 to 2 simulation sound 
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crossing times) pressure and density histograms for two fully turbulent simulations with 
driving wavenumber fcf or = 2 and rms Mach number M — 1, at resolutions of 512 2 and 1024 2 
(fig. 3). We see that the right sides (high values) of the histograms are perfectly converged 
at 512 2 , while the left sides (low values) are approximately so, with the same relative number 
of points as the 1024 2 run in density and pressure intervals within a factor of 3 from the 
histogram maximum, and deviations by factors no larger than ~ 3 in more distant density or 
pressure intervals. From all of the above results, we therefore adopt N = 512 as a compromise 
between acceptable resolution, and the ability to perform the numerous simulations needed 
for this study. 

Note that with a box of 100 pc on each side and N = 512, the smallest resolved scale, 
even neglecting the effects of numerical diffusion, is ~ 0.2 pc. This is larger than the typical 
size (~ 0.1 pc) of the cold structures generated by TI and thus our simulations probably 
overestimate the sizes of those cloudlets that are formed by the instability rather than by 
larger-scale, coherent turbulent compressions, and certainly do not resolve their internal 
structure. However, our interest here is in statistical quantities, such as the distribution of the 
pressure values at every density interval, and the fraction of the volume occupied by gas with 
a given pressure. The fact that the high-value sides of the density and pressure histograms 
at 512 2 and 1024 2 resolutions are nearly identical suggests that this information deos not 
require resolving the tiniest structures, and is accurately captured by our simulations. 



4. Results 

We now describe the results of 10 simulations aimed at characterizing the effect of 
velocity fluctuations on the thermal pressure distribution in a thermally bistable flow. The 
simulations analyzed in the next two sections are performed in 2D, with a resolution of 
512 grid points per dimension, whereas the simulations described in §5.2 are performed in 
two and three dimensions, with 256 grid points in each direction. The forcing is applied at 
wavenumbers of k ioi = 2, 4, 8 and 16, implying a driving scale of 50, 25, 12.5 and 6.25 pc 
respectively, and has the necessary amplitude to induce turbulent motions with rms Mach 
numbers M ~ 0.5, 1.0 and 1.25. This is intended to represent realistic values of the Mach 
number in the warm phase. Of course, actual local Mach numbers can be much higher, as 
the temperature is generally lower in higher density gas. 

Initially, all our simulations are at rest and have a uniform density (no = 1 cm -3 ) and 
temperature (To = 2400K). In order to allow the simulations to reach a stationary regime, 
we evolve the simulations for at least two crossing times across the turbulence driving scale, 
given by t /Mk ior . Specifically, cases with k {or > 4 are evolved for two code time units (2t ), 
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while those with kf OT = 2 are evolved for 4t - These times represent several cooling times at 
T , given by 

where cy is the specific heat at constant volume and p eq is the unstable equilibrium density. 



4.1. The effect of the driving scale 

We first discuss the response of the pressure to changes in the driving wavenumber fcf or , 
and so in this section we restrict ourselves to rms Mach numbers M pa 1. Recalling that M 
is not directly an input parameter of the simulations, but a result of the energy input rate, 
the actual value of M differs slightly from the target value. Specifically, its actual values are 
M = 0.97, 0.96, 0.92 and 0.95 for the runs with fcf or = 2, 4, 8 and 16, respectively. 

Figure 4 shows two-dimensional histograms for these runs, giving the number of grid 
points in the simulation in a given (P, n) bin. The contours are logarithmic and are set at 
10%, 30%, 50%, 70% and 90% of the logio of the maximum value of the two-dimensional 
histogram for each simulation. The histograms are computed at t — 1.5£o (resp. t = 3.0to) 
for simulations with kf OT > 2 (resp. fcf or = 2). It is clearly seen that as the driving scale 
decreases (fcf or increases), the distribution of points shifts away from the thermal equilibrium 
curve (denoted by the broken solid line), and towards adiabatic behavior (with slope 5/3). 
It is also seen that at low densities a substantial fraction of the points lies below the thermal 
equilibrium curve. That is, they have probably been cooled by negative P dV work, and 
have not had time yet to warm back up by the background heating. Finally, an interesting 
branch of points seems to lie on the extension of the equilibrium curve for the dense gas, but 
at densities corresponding to the diffuse gas. 

These trends can also be seen in figure 5, which shows temporally- averaged pressure 
histograms computed in three density ranges n c j ' \/2 < n c < V2n c , with n c = 0.1, 1.0 and 
10.0 cm~ 3 , corresponding to the warm, unstable and cold ranges. It can be seen that for low 
densities the most probable pressure P(iV max ) is in general lower than the thermal equilib- 
rium pressure at n c (P eq , denoted in the figure by the vertical lines), and shifts progressively 
farther away from it as fcf or increases. At densities corresponding to the thermally unsta- 
ble range, the most probable pressure increases with increasing kf OT , and passes from being 
smaller to being larger than P eq . Finally, at n c = 10.0 cm -3 , the four histograms peak close 
to P cq , although the height decreases rapidly with increasing kf OT . Moreover, a high-pressure 
tail is present, which becomes more populated and more extended as k{ or is increased. Fi- 
nally, we note that, in all histograms displayed in figure 5, the dynamic range is larger than 
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an order of magnitude. 

The time-averaged pressure and density histograms for the whole field of all the sim- 
ulations in this group are shown in figure 6. The pressure histograms span 3-4 orders of 
magnitude, with both their width and P(iV max ) increasing as kf OT increases. Also, the his- 
tograms become more skewed, lifting their low-P side. The density histograms are bimodal, 
but as fcf or increases, the bimodality becomes less pronounced and the histogram becomes 
narrower. 



4.2. The effect of the Mach number 

We now turn to the effect of the rms Mach number on the pressure distribution. To 
this effect, in this section we present results from two sets of three simulations each, with 
M ~ 0.5, 1.0 and 1.25, and kf or = 2 and 16. The actual values of M in these simulations are 
M = 0.50, 0.97 and 1.25 at k ioi = 2, and M = 0.56, 0.95 and 1.25 at k ior = 16. 

The two-dimensional histograms in the pressure-density space are shown in figures 7 
and 8 for kf OT = 2 and 16, respectively. The contours are set at the same levels as in §4.1. 
The dynamic range of both density and pressure is seen to increase with increasing M for 
both driving scales. For kf or = 2, a slight variation of the mean slope of the distribution of 
points can be easily seen, while a clear steepening of the mean slope is observed at kf ov = 16. 
This variation is summarized in fig. 13 for all runs. Finally, the tendency of the pressure 
distribution to shift away from thermal equilibrium for smaller driving scales, reported in 
the previous section, is seen for all three values of M. 

The pressure histograms for specific density intervals show a variety of behaviors. For 
/cf or = 2 (fig. 9), the histogram centered on n c = 0.1 cm -3 shifts from being narrow and 
peaking very close to P cq at M = 0.5 to a wider distribution peaking below P cq at M — 1 
and 1.25. The histograms for the unstable gas at n c = 1 exhibit almost no shift of their peak 
P(N max ), which is located a half order of magnitude below P eq , as M is increased, although 
the high-P tail becomes higher, with a shallower slope, and extends to higher pressures. In 
the high-n range, the histograms show almost no change in P(iV max ), but the high-P tail 
becomes progressively shallower, and extends up to higher values, as M is increased. 

For k ior = 16 (fig. 10), the histograms at n c = 0.1 cm -3 and n c = 1 cm -3 show little 
variation in shape and in P(iV max ), although the histogram at n c — 0.1 almost doubles its 
height, indicating that a higher fraction of the volume in the simulation is occupied by diffuse 
gas as M is increased. This is compensated by a decrease in the volume occupied by the 
dense, n c = 10 cm -3 gas, which moreover experiences an increment of roughly one and a 
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half orders of magnitude in its typical pressure P(A" max ). 

In the pressure histograms for the whole simulations (fig. 11 left) at k {or = 2, the 
different behavior between the run with M ~ 0.5 and those with M ~ 1.0 and M ~ 1.25 is 
also evident. In the former case the histogram is narrower and close to lognormal, although 
with a high-P tail for P > 10 4 K cm -3 . As M increases, P(N ma , x ) shifts to slightly lower 
pressures, and the high-P tail reaches closer to the histogram peak, lifting the entire high-P 
side of the histogram and causing it to approximate a power law. On the other hand, the low- 
P side of the histogram widens with increasing M, but never seems to lose its approximately 
lognormal shape. 

In the density histograms for kf OT = 2 (fig. 11 right), the distinction between M ~ 0.5 
and the other two cases is also noticeable. The bimodal shape of the histogram becomes less 
pronounced as M is increased, with n(N max ) shifting slightly towards lower densities. The 
histogram width increases with increasing M, at least in the range explored. 

For fcf or = 16 (fig. 12, left), the high-P branch in the pressure histogram decays faster 
than the low pressure one at all values of M, and its slope is almost independent of M, 
although P(A max ) shifts to higher pressures. In this case again the histogram becomes 
broader with increasing M, but mainly because the low-P tail is lifted and extends to lower 
pressures. 

Concerning the density histograms, similarly to the case with fcf or = 2, for kf OT = 16, 
an increase in M leads to a less bimodal density histogram (fig. 12, right), but in this case, 
the histogram broadens with increasing M, although mainly by lifting its low-n tail. A 
comparison between the density histograms shown in figures 11 and 12 confirms the fact, 
discussed in previous section, that the density distribution becomes narrower with increasing 
hoi- 

5. Discussion 

5.1. A unified physical picture 

The main results of §4 can be summarized as follows: (1) The distribution of points 
in the P-n diagram widens and steepens as either M or fcf or are increased. (2) The mean 
pressure in a given density interval drifts away from the equilibrium value P eq as fcf or is 
increased, moving towards P > P eq for the dense gas, and towards P < P eq for the diffuse 
gas. (3) The pressure histograms in these density ranges as well as the global pressure 
histograms are generally skewed, and tend to reverse their skewness with increasing kf OT . (4) 
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The pressure histograms in specific density ranges increase their width as M is increased. 
(5) The global pressure histograms develop one near-power-law side and one near-lognormal 
side. For low fcf or , the near-power-law develops on the high-P side, while for high fcf or the 
near-power-law side develops at low P. In general, the slope of the power law flattens as M 
is increased. 

Most of these results can be understood simply as a consequence of the turbulent cross- 
ing time becoming shorter in local compressions as either the Mach number M or the driving 
wavenumber kf OT are increased, creating a larger fraction of compressions that evolve closer 
to adiabatically, and thus temporarily drift away from thermal equilibrium. Indeed, for our 
choice of parameters, the cooling and sound crossing times in the unstable gas are equal at a 
scale A cq ~ 4 pc. For the warm medium in thermal equilibrium, this scale increases to ~ 23 
pc. The scale A eq also applies for equality of the turbulent crossing time and the cooling time 
for Mach-1 motions. Below this scale, classical isobaric perturbations would evolve nearly 
isobarically, because condensation occurs on roughly the cooling time, which is longer than 
the time needed to restore pressure balance (the sound crossing time). However, velocity 
perturbations below this scale generate perturbations that approach adiabatic behavior as 
their Mach number increases, because in this case the externally-applied turbulent compres- 
sion exerts PdV work on the fluid parcel, heating it on timescales shorter than the cooling 
time. Adiabatic perturbations are stable to first order (Field 1965), and imply that the 
pressure increases with increasing density. The net behavior in transonic flows is expected 
to be intermediate between thermal equilibrium and adiabaticity. 

The tendency towards adiabatic behavior causes a progressive increase in the slope (or 
effective polytropic exponent, 7 C ) of the ensemble of points in the P-n diagram (fig. 13, 
left panel). This causes a decrease in the mean pressure of the diffuse gas, and an increase 
in the mean pressure of the dense gas, because the point distribution is centered in the 
unstable range. For the dense gas, this furthermore causes a tendency towards producing 
flatter-topped histograms, because the fraction of high-pressure zones increases, but the 
short cooling time at those densities always produces a significant fraction (most frequently 
a majority) of points near P cq . Only for the case with the highest M and fcf or does the 
peak of the pressure histogram for the dense gas shift to the high-P part of the histogram. 
Nevertheless, the fraction of zones in the simulations with P > 10 4 K cm -3 is seen to increase 
monotonically with either M or kf OT (fig. 13, right panel). The fraction of the total number of 
cells with n > 7.1 cm" 3 and P > 10 4 K cm 3 increases from 0.07% at kf OT = 2 and M = 0.5 
to 69% at k ioi = 16 and M = 1.25. 

In addition, larger rms Mach numbers are known to cause wider density PDFs (Padoan 
et al. 1997; Passot & Vazquez-Semadeni 1998), with the amplitude of the density fluctuations 
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depending on 7 C (Vazquez- Semadeni et al. 1996). As a consequence, wider pressure PDFs 
are also expected for stiffer-than-isobaric (i.e., 7 C > 0) behavior, in which the pressure is 
positively correlated with the density. 

For the density PDF, Passot & Vazquez- Semadeni (1998; see also Scalo et al. 1998; 
Nordlund & Padoan 1999) showed that the density PDF is lognormal for isothermal flows 
(7e = 1; see also Vazquez-Semadeni 1994), but develops a power-law tail at high densities 
for 7c < 1, and at low densities for 7 e > 1. A similar trend is also observed here in response 
to the resulting effective polytropic exponent: simulations with /cf or = 2 have 7 C < 0.5 (fig. 
13, left), and their resulting density PDFs are skewed to the left, with shallower high-n tails 
that extend further from the PDF peak than the low-n side (fig. 11 right), although with 
signatures of the bimodality associated with the thermal bistability for the cases with the 
lowest-M. Instead, simulations with k {or = 16, have density PDFs in which the low-n side is 
shallower and generally more extended, although again with signatures of bimodality on the 
high-n side at low M (fig. 12 right). Thus, the hf OT = 2 runs generally behave as if having 
7e < 1, while k{ or = 16 runs behave as if having 7 e > 1. 

Interestingly, this dependence of the density PDF on 7 C is also apparent in the pressure 
histograms, and in fact it is even more pronounced, as can be seen in figs. 11 and 12 (left 
panels). This is somewhat surprising, given the large scatter of the P-n points in any given 
density interval. Naively, one would expect that such a large scatter would preclude any 
copying of the density-PDF features into the pressure histogram. Indeed, most features of 
the distribution do not survive the change of variable from density to pressure. This is the 
case of the scaling of the histogram width or the shift in position of the histogram peak with 
Mach number. Nevertheless, the development of a near power-law tail depending on the 
value of 7 e does seem to be preserved, and even amplified, in the pressure histogram. 

5.2. Three-dimensional tests 

The results from §§4.1 and 4.2 have been obtained in two-dimensional (2D) simulations 
exclusively, and it is thus important to determine whether these results are expected to persist 
in three dimensions (3D). The distinction between the 2D and 3D cases has been extensively 
discussed by Vazquez-Semadeni (1994) and Avila- Reese & Vazquez-Semadeni (2001). These 
authors have argued that the distinction is less pronounced in the highly compressible case, 
in which shocks are an important ingredient in the dynamics. This is because shocks are 
essentially one-dimensional structures, independently of the dimensionality of the global flow. 

In order to confirm this expectation in the case of our particular problem, in this section 
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we present a comparison between selected cases in 2D and 3D, using simulations with 256 
grid points in each direction. The 3D simulations are performed using a parallel version of 
the code. We first consider two simulations with fcf or = 2 and M ~ 0.9, one in 2D and the 
other in 3D. The precise Mach numbers are 0.90 and 0.85, respectively. Figure 14 shows 
the total histograms of pressure (left) and density (right). There we see that the pressure 
distribution of the 3D run is systematically shifted to higher values, by a factor < 2, with 
respect to that of the 2D one. On the other hand, the density distributions of the two runs 
coincide at high densities, but the 3D case has its peak and its low-density side again shifted 
to higher densities by about a factor of 2. 

We speculate that these effects may be due to the fact that the density peaks are built 
by compressions of varying dimensionality, up to the dimensionality of the flow. The higher 
fraction of low-density regions in 2D can be understood as a consequence that a peak of 
a given average density and size contains a larger fraction of the total mass in 2D than in 
3D. Thus, the "voids" surrounding the peaks are more heavily evacuated, giving a higher 
fraction of low-density zones and histograms that extend to lower values of the density. This 
typically lower density of the voids in 2D can also explain the typically lower values of the 
pressure at the low-pressure side of the distribution. On the other hand, for a density peak 
formed by a compression at a given characteristic velocity R = dR/dt, the density varies 
more rapidly in 3D than in 2D. That is, assume that p oc M.Rr m ', where M. and R are 
respectively the mass and radius of the compressed parcel, and m = 2 in 2D and 3 in 3D. 
Then p = —mM.R~ m ~ l R, and the density rate of change is larger in 3D than in 2D at a 
given compression velocity R. Thus, in 3D the characteristic time for variation of the density 
is comparatively shorter than the cooling time in 3D, and the behavior should be slightly 
closer to adiabatic, explaining the higher transient pressures in 3D at high densities. 

Nevertheless, the differences in the histograms in 2D and 3D are relatively minor, prob- 
ably because the occurrence of high-dimensional compressions must be a relatively rare event 
in comparison with one-dimensional ones (shocks). The pressure distribution in 3D is shifted 
by factors not larger than 2 in the 3D case, and in fact the P-n relation is very similar in the 
2D and 3D cases, as seen from two-dimensional histograms in the pressure-density space (fig. 
15), with the main difference being that the distribution of points in the 3D case extends 
slightly farther above the thermal equilibrium curve at densities n ~ 1 cm" 3 . The least 
squares slopes of the distributions are also very similar, at 0.35 and 0.38 for 2D and 3D, 
respectively. 

The trends described in previous sections for the pressure and density distributions 
resulting from 2D simulations as M and kf or are varied are also maintained in the 3D case. 
In fig. 16a and b we respectively show the pressure and density histograms resulting from 
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three simulations with (M, kf OT ) = (0.85,2), (0.90,8), and (1.35,2). It can be observed that 
as fcf or increases, the width of the pressure distribution and the value of P(iV max ) increase 
while the density histogram becomes narrower and with a less pronounced bimodality. On 
the other hand, when the value of M is increased, the pressure distribution widens noticeably, 
while the density distribution widens marginally; also, the bimodality of the latter is slightly 
less pronounced. 

Finally, the fraction of cells with n > 7.1 cm~ 3 and P > 10 4 K cm -3 for the 3D 
simulation with fcf or = 2 and M = 0.85 is 1.3%, while that obtained in the 2D run with 
fcf or = 2 and M = 0.97 is 1.2%, again showing a high consistency between the 2D and 3D 
cases. 

We conclude that the results obtained from the 2D simulations presented in §§4.1 and 
4.2, as well as the discussion from §5.1, still hold in 3D. 

5.3. Relation to previous work 

Our results support the scenario that in a turbulent, thermally bistable flow, there 
exists a fraction of fluid parcels that are out of thermal equilibrium, and that this fraction 
depends on the local ratio r\ of the turbulent crossing time to the cooling time. Recently, 
Wolfire et al. (2003) have estimated this ratio for the warm neutral medium (WNM), using 
an approximation for the turbulent crossing time given by t s h oc k ~ A s /c s (i.e., the "mean 
time between shocks"), where c s is the sound speed and A s is the scale at which the typical 
turbulent velocity difference equals c s . We refer to A s as the "sonic" scale. Wolfire et al. 
(2003) estimated A s ~ 200 pc for the WNM, and ~ 0.3 pc for the cold neutral medium 
(CNM), although they warned that this is a very uncertain quantity. With these estimates, 
they found 77 ~ 0.3-0.9 for the WNM, concluding that non-equilibrium temperatures should 
often be found in the WNM, in agreement with observations (Dickey et al. 1977; Kalberla 
et al. 1985; Spitzer & Fitzpatrick 1995; Fitzpatrick & Spitzer 1997; Heiles 2001; Heiles & 
Troland 2003; Kanekar et al. 2003) and previous numerical studies (Gazol et al. 2001). 

This conclusion is also consistent with our present results. Adopting their value of 
~ 200 pc for A s , and a Kolmogorov velocity dispersion scaling law Av oc A 1 / 3 , appropriate 
for incompressible turbulence, we see that the rms Mach number with respect to the WNM 
at a scale of our simulation boxes (100 pc), should be ~ 2 -1 / 3 0.8. Thus, based on 
their estimates, the WNM on scales of 100 pc should be bracketed by our simulations with 
M = 0.5 and M — 1, at the largest driving scales (hf OI = 2), with the M — 1 case being the 
most relevant. 
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However, it is possible that the above regime, arrived at through the considerations of 
Wolfire et al. (2003), still somewhat underestimates the role of turbulence. Those authors 
noted that, because of the scaling of the turbulent velocity with size, below the sonic scale 
the medium should roughly be in pressure equilibrium. However, this does not guarantee 
that thermal instability will be fully unimpeded in the generation of the density structures 
below this scale. As discussed by Vazquez- Semadeni et al. (2003), in a turbulent, thermally- 
bistable flow, there are three competing timescales: the sound crossing time r s , the turbulent 
crossing time r t , and the cooling time r c . In the standard linear analysis (Field 1965), the 
largest growth rate of the instability is given by the inverse of the cooling time, which is 
independent of scale in the linear regime. This largest growth rate occurs at an intermediate- 
wavelength regime Af <C A < A eq , where Af is the Field length and A eq is the scale at which 
r s ~ T c ( c -f- §5.1). On the other hand, the sound and turbulent crossing times do depend 
on scale, with r s oc A and r t oc (assuming an incompressible Kolmogorov spectrum), 
or r t oc A 1 / 2 (assuming a highly compressible, Burgers-like spectrum, appropriate for the 
unstable range). So, even though indeed the ratio of sound-to-turbulent crossing times 
becomes progressively smaller with decreasing scale size, the ratio rj = r t /r c also becomes 
smaller. This suggests that nearly incompressible, shearing turbulent fluctuations may have 
time to disrupt Tl-induced condensations before they grow, at least partially. This effect 
may be enhanced in the presence of magnetic fields, which increase the solenoidal (shearing) 
fraction of the turbulent kinetic energy (Vazquez- Semadeni et al. 1996). Indeed, even in our 
weakest-turbulence (largest-77) simulations (M = 0.5, fcf or = 2), the bimodality of the density 
PDF (caused by the thermal bistability) is moderate (fig. 11 right panel), and the global 
effective polytropic exponent is already positive across the thermally unstable range (fig. 13, 
left panel)). 

Our results also have the implication that the development of a high-P power-law tail 
is not exclusive to supernova-driven models, such as the McKee & Ostriker (1977) model, 
and that pure advective driving can generate approximate power-law tails in low-7 e flows at 
sufficiently high Mach numbers. On the other hand, noting that we also obtain lognormal- 
like shapes in the high-P sides of the PDFs for large fcf or suggests that the formation (or 
not) of power-law tails in the pressure PDF depends on j e . 

Finally, the fact that the pressure histogram partially copies features of the density 
histogram, in particular the production of power-law tails as a function of j e , gives validation 
to the approach used by Mac Low et al. (2004), who applied a change of variable from density 
to pressure in order to understand the functional form of the pressure PDF from the known 
shape of the density PDF. On first account, such a procedure might be questioned because 
of the large scatter exhibited by the pressure at any given density, but our results suggest 
that it may be valid on average. 
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Observationally, Jenkins (2004) has recently presented histograms of the pressure in 
observations of CI, which can in principle be compared with our results, as well as with 
pressure histograms from full ISM simulations, like those of Mac Low et al. (2004) and 
de Avillez & Breitschwerdt (2004). However, at this point such a comparison is ambiguous, 
because Jenkins (2004) had to assume an effective polytropic exponent in order to correct for 
an over-representation of dense points in his observational sample. Since our results indicate 
that 7 e determines the shape of the pressure histogram, it appears difficult to disentangle 
this physical role of 7 e from its possible role in biasing the observed pressure distribution. 
A detailed attempt to perform such a comparison will be presented elsewhere, but here we 
just note that our fraction of points with P > 10 4 K cm -3 for the case M = 0.5, k ior = 2, of 
order 0.07%, appears consistent with the fraction reported by Jenkins (2004) for approximate 
thermal equilibrium in the cold medium, also ~ 0.1%. 



6. Conclusions 

In this paper we have carried out a systematic study of the effect of turbulent velocity 
fluctuations, characterized by their rms Mach number M and the energy injection wavenum- 
ber, fcf or , on the thermal pressure distribution in a thermally bistable flow. To this end we 
have performed a large number of 2D simulations varying those two parameters in 100-pc 
boxes with random turbulent driving generated in Fourier space, which allows precise control 
of the parameters. A few test cases in 3D suggest that the 2D results are valid in 3D as well. 
Our results are consistent with the picture that as either of these parameters is increased, 
the ratio of turbulent crossing time to cooling time decreases, causing a departure from 
thermal equilibrium, and an approach towards an adiabatic behavior. This translates into 
an increase of the effective polytropic index j e , as measured by the slope of the distribution 
of points in the pressure- density diagram, in turn creating a population of underpressured 
zones in the diffuse gas, and of overpressured zones in the dense gas, with respect to the 
thermal-equilibrium value of the pressure, P eq . In particular, the fraction of zones with den- 
sities n > 7.1 cm" 3 and with P > 10 4 K cm~ 3 increases from roughly 0.1% at fcf or = 2 (a 
driving scale of 50 pc) and M = 0.5 to roughly 70% for fcf or = 16 (a driving scale of 6.25 pc) 
and M = 1.25. In particular, for M — 1 and k lor = 2, this fraction is ~ 1%, similar to the 
value reported by Jenkins (2004) from a survey of the fine-structure excitation of CI on the 
Glactic plane. 

In all cases, the dynamic range of the pressure in any given density interval is larger 
than one order of magnitude, and the total dynamic range, summed over the entire density 
range, typically spans 3-4 orders of magnitude. The total pressure histogram widens as the 
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Mach number is increased, and moreover develops near-power- law tails at high (resp. low) 
pressures when 7 e < 0.5 (resp. 7 C > 1), which occurs at k ior = 2 (resp. k {oT = 16) in our 
simulations. The opposite side of the pressure histogram decays rapidly, in an approximately 
lognormal form. The development of power-law tails in the pressure PDF is analogous to, 
and in fact more pronounced than, that observed in the density PDF, suggesting that the 
average value of the pressure is sufficiently correlated with the density as to reflect the same 
dependence of its histogram with j e , in spite of the large scatter of the pressure at any given 
density. This may validate approaches attempting to obtain the pressure PDF from that of 
the density via a change of variables. 
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Fig. 1. — Temporal growth of sinusoidal density perturbations with an initial amplitude of 
2.5% and a period of 50 pc (upper left), 25 pc (upper right) 12.5 pc (lower left) and 6.25 pc 
(lower right). The dashed, dotted and solid lines correspond to resolution of 256 2 , 512 2 and 
1024 2 , respectively. The thick straight lines indicate the slope of the theoretical growth rate 
at the corresponding scale. 
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Fig. 2. — Temporal growth of sinusoidal density perturbations with an initial amplitude of 
2.5% and a period of 6.25 pc (left) and 12.5 pc (right) for simulations including an initial 
velocity perturbation at amplitude equivalent to Mach number M — 1.0 with respect to the 
unstable medium at T = 2400 K. The dashed, dotted and solid lines correspond to resolutions 
of 256 2 , 512 2 and 1024 2 , respectively. 
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Fig. 3. — Comparison of the total pressure (left) and density (right) histograms for sim- 
ulations with resolutions 512 2 (dotted lines) and 1024 2 (solid lines). The histograms are 
normalized to the total number of points, and averaged over the time interval 1.1 < t < 2, 
where the time unit is the simulation sound crossing time at T = 10 4 K 
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Fig. 4. — Thermal pressure-density relation for simulations with M ~ 1 and fcf or = 2 (upper 
left), fcf or = 4 (upper right), fcf or = 8 (lower left) and fcf or = 16 (lower right). The solid line in 
each panel denotes the thermal-equilibrium pressure. 
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Fig. 5. — Temporally-averaged pressure histograms for simulations with M ~ 1 and kf ov = 2 
(solid line), fcf or = 4 (dotted line), kf or = 8 (dashed line) and fcf or = 16 (dashed- dotted line). 
The time averaging is as in fig. 3. The histograms are computed over logarithmic density 
intervals n c / y/2 < n c < V~2n c , centered at n c — 0.1 cm~ 3 (upper left panel), n c = 1.0 cm~ 3 
(upper right panel), and n c = 10.0 cm~ 3 (lower left panel). The vertical lines denote the 
thermal-equilibrium pressure P eq at n c . Lower right panel: most probable pressure P(N max ) 
for the histograms centered on n c = 0.1 cm" 3 (plus signs) and n c = 1.0 cm" 3 (asterisks) . The 
horizontal lines denote the values of P eq (n c = 0.1 cm" 3 ) (solid line) and P eq (n c = 1 cm -3 ) 
(dotted line). 
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Fig. 6. — Total pressure (left) and density (right) histograms for simulations with M ~ 1 and 
fcfor = 2 (solid line), kf or = 4 (dotted line), kf OI = 8 (dashed line) and fcf or = 16 (dashed- dotted 
line). The histograms are time- averaged as in fig. 3. 
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Fig. 7. — Thermal pressure-density relation for simulations with /c for = 2 and M ~ 0.5 (upper 
panel), M ~ 1 (center panel) and M ~ 1.25 (lower panel) . The solid line in each panel shows 
the thermal-equilibrium pressure. 
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Fig. 8. — Thermal pressure-density relation for simulations with k iov = 16 and M ~ 0.5 
(upper panel), M ~ 1 (center panel) and M ~ 1.25 (lower panel). The solid line in each 
panel shows the thermal-equilibrium pressure. 
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Fig. 9. — Pressure histograms over specific density intervals as in fig. 5, but for simulations 
with k {oT = 2 and M ~ 0.5 (solid line), M ~ 1.0 (dotted line), and M ~ 1.25 (dashed line). 
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Fig. 10. — Pressure histograms over specific density intervals as in fig. 5, but for simulations 
with fcf or = 16 and M ~ 0.5 (solid line), M ~ 1.0 (dotted line). 
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Fig. 11. — Total time-averaged pressure (left) and density (right) histograms as in fig. 6 but 
for simulations with kf OT = 2 at M ~ 0.5 (solid line), M ~ 1.0 (dotted line), and M ~ 1.25 
(dashed line). 
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Fig. 12. — Temporally-averaged pressure (left) and density (right) histograms as in fig. 6 
but for simulations with kf OT = 16 and M ~ 0.5 (solid line), M ~ 1.0 (dotted line), and 
M ~ 1.25 (dashed line). 
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Fig. 13. — Left: Least squares slope of the distributions of points in figs. 4, 7 and 8. The solid 
line shows the variation with the driving wavenumber kf OT , indicated by the lower horizontal 
axis, at fixed rms Mach numbers M — 1. The dotted and dashed lines show the variation 
with M, indicated by the upper horizontal axis, at a given kf or , indicated by the label next 
to each line. Right: Fraction of grid cells with P > 10 4 K cm" 3 and for the n > 7.1 cm -3 in 
figs. The line coding is as in the left panel. 
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Fig. 14. — Pressure (left) and density (right) histograms for simulations with k{ OT = 2 and 
M = 0.90 in 2D (solid line), and M = 0.85 in 3D (dotted line). The histograms are 
normalized to the total number of points and for the 2D simulations they are averaged over 
the time interval 3.2 < t/t < 4. 
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Fig. 15. — Thermal pressure-density relation for three-dimensional simulations with fcf or = 2 
and M = 0.90 in 2D (left panel) and M = 0.85 in 3D (right panel). The solid line in each 
panel shows the thermal-equilibrium pressure. 
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Fig. 16. — Total pressure (left) and density (right) histograms for 3D simulations with 
M = 0.85 and kf or = 2 (solid line), M = 0.85 and k loi = 8 (dotted line), and M = 1.35 and 
fcf or = 2 (dashed line). The histograms are normalized to the total number of points. 



